Beyond Mean-Field Theory for Attractive Bosons 
under Transverse Harmonic Confinement 
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We study a dilute gas of attractive bosons confined in a harmonic cylinder, i.e. under cylindric 
confinement due to a transverse harmonic potential. We introduce a many-body wave function 
which extends the Bethe ansatz proposed by McGuire (J. Math. Phys. 5, 622 (1964)) by including 
a variational transverse Gaussian shape. We investigate the ground state properties of the system 
comparing them with the ones of the one-dimensional (ID) attractive Bose gas. We find that the 
gas becomes ultra ID as a consequence of the attractive interaction: the transverse width of the 
Bose gas reduces by increasing the number of particles up to a critical width below which there is 
the collapse of the cloud. In addition, we derive a simple analytical expression for the simmetry- 
breaking solitonic density profile of the ground-state, which generalize the one deduced by Calogero 
and Degasperis (Phys. Rev. A 11, 265 (1975)). This bright-soliton analytical solution shows near 
the collapse small deviations with respect to the 3D mean-field numerical solution. Finally, we show 
that our variational Gauss-McGuire theory is always more accurate than the McGuire theory. In 
addition, we prove that for small numbers of particles the Gauss-McGuire theory is more reliable 
than the mean-field theory described by the 3D Gross-Pitaevskii equation. 



PACS Numbers: 03.75.Kk 



I. INTRODUCTION 



A three-dimensional (3D) dilute Bose gas is well de- 
scribed by the mean-field Gross-Pitaevskii theory [1,2], 
while in the ID regime quantum fluctuations become 
important and an appropriate theoretical treatment re- 
quires beyond mean-field approaches [3]. Recently, we 
have investigated, for a bosonic cloud of atoms, the 
crossover from 3D to ID induced by a strong harmonic 
confinement in the transverse cylindric radial direction 
[4-7]. In particular, by using a generalized Lieb-Liniger 
theory [8,7], based on a variational treatment of the 
transverse width of the repulsive Bose gas, we have ana- 
lyzed the transition from a 3D Bose-Einstein condensate 
to a ID Tonks- Girardeau gas of impenetrable bosons. 

In this paper we consider the case of an attractive Bose 
gas under cylindric transverse harmonic confinement. We 
introduce a trial wave fuction that is a variational ex- 
tension of the ID exact Bethe ansatz [10] proposed by 
McGuire [11] for the ID gas of bosons with attractive 
contact interaction. We show that, contrary to the ID 
theory where the transverse width is constant and equal 
to the characteristic length of harmonic confinement, our 
theory predicts that the transverse width of the gas re- 
duces by increasing the interatomic strength up to a criti- 
cal value for which there is the collapse of the system. We 
investigate also the solitonic axial density profile of the 
Bose gas with a fixed center of mass, comparing it to the 
ID profile and also to the 3D mean-field Gross-Pitaevskii 
theory. 



II. TRANSVERSE GAUSSIAN ANSATZ 

The Hamiltonian of a gas of N interacting identical 
Bose atoms confined in the transverse cylindric radial 
direction with a harmonic potential of frequency uj±_ is 
given by 
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E nr„r,), (1) 



where V{ri,rj) is the inter-atomic potential. In the 
Hamiltonian we use scaled units: energies are in imits 
of the energy huj± of the transverse confinement and 
lengths in units of the characteristic harmonic length 
ax = {h/{mu;^)y/^. 

The determination of the A^-body wave function 
^'(ri, ...,rjv) that minimizes the energy 



E = {■^\H\'i') = / d^ri...(frN 



(2) 



of the system is a very difficult task. Nevertheless, due 
to the symmetry of the problem, a variational trial wave 
function can be written in the form 



«'(ri,...,rAr) = /(zi,. 



N 
i=l 



exp 
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where cr is a variational parameter of the Gaussian trans- 
verse wave function [13], that gives the effective trans- 
verse length of the Bose gas. In addition, by considering 
a dilute gas with a mean particle spacing much larger 
than the interaction radius we set 



T/(r„r,)=r<5(3)(r,,_r^.) 



(4) 
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where T = Anas/a^, with as the s-wave scattering length 
of the inter-atomic potential. This pseudo-potential gives 
the correct dilute gas limit and a well-posed variational 
problem by choosing a smooth trial wave function. Under 
transverse harmonic trapping a confinement-induced res- 
onance has been predicted by Olshanii [9] at \as\/a± ^ 1, 
i.e. when the absolute value of the highest bound-state 
energy of a realistic inter-atomic potential approaches 
the confining transverse energy. Therefore the pseudo- 
potential of Eq. (4) can be used in the range \as\/a± <gC 1 , 
a regime where the effects of confinement-induced reso- 
nance are negligible. 

By inserting Eq. (3) into Eq. (2), using Eq. (4) and 
integrating over x and y, the total energy reads 



E = Ez + E± , 
where the longitudinal axial energy is 
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It is important to stress that the ID Hamiltonian that 
appears in the previous expresssion is exactly the Hamil- 
tonian studied by Lieb and Liniger [8] for positive inter- 
action strength (r/(27r(7^) > 0) and by Mc Giurc [11] and 
by Calogero and De Gasperis [12] for negative interaction 
strength (r/(27ra-^) < 0). The transverse radial energy is 
instead given by 



(7) 



As previously stressed, with the ansatz of Eq. (3) we have 
recently analyzed [6] the repulsive case (a^ > 0) by using 
the Licb-Liniger exact result [8] for the longitudinal axial 
energy. Here we consider the attractive case (a^ < 0) and 
set r/27rcr^ = -7/0-^ with 7 = 2\as\/a±. 



III. GAUSS-MCGUIRE ANSATZ 

McGuire [11] proposed the following Bethe ansatz [10] 
for the axial many-body wave function 

f{zi,...ZN) =Cn n expf-i^l^i-^r^l j , (8) 

l<i<j<N \ ^ / 

where Cjv is the normalization constant and ^jo^ is the 
strength of the contact (5 interaction in the ID Hamilto- 
nian of Eq. (6). According to the exact result of Mcguire 
[11] the axial energy of Eq. (6) reads 



1 -v^ 
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Using this expression the total energy E of our 3D system 
can be rewritten as 



E: 



■l^(iV2-l)7V+i(^+aV. (10) 



In the weak-coupling ID limit, where the transverse 
width (T ~ 1, one recovers the ID result obtained by 
McGuire [11] plus the constant transverse energy, that is 
1 in units of Tiwj,. 

In our approach the energy depends on the variational 
parameter a. The minimization of the energy E with 
respect to a gives the equation 



(11) 



One easily finds that this algebric equation admits real 
solutions with dPE/da'^ > if and only if 



1 

3174 



< fj < 1 . 



(12) 



Below the value a = 1/3^/^ there are no stable solutions 
and the ground state becomes the collapsed state, i.e. 
the configuration with cr = and energy E = —00. The 
critical strength, corresponding to cr = 1/3^^*, is 



(13) 
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FIG. 1. Comparison between our variational Gauss- 
McGuire theory (GMG) and the McGuire theory (MG) 
for the attractive Bosc gas in a harmonic cylinder. Up- 
per panel: transverse width a. The unstable branch 
{(PE/da"^ < 0) of Eq. (11) is shown as a dotted line. 
Lower panel: energy per particle E/N. 7 = 2|as|/o_L is 
the inter-atomic strength and N is the number of parti- 
cles. 

In the upper panel of Fig. 1 we plot the trans- 
verse width a obtained from Eq. (11) as a function of 
7(Ar2 _ ;^^i/2 ^gQjjjj jjj^g^ While in the ID case a is 
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constant (dashed line), our variational Gauss-McGuire 
(GMG) method shows that the attractive Bosc gas is 
actually ultra ID: the width a decreases by increasing 
the interaction strength up to the collapse. In the lower 
panel of Fig. 1 we plot the energy per particle E/N as 
a function of ^{N"^ — 1)^/^. We compare our variational 
energy given by Eqs. (10,11) with the McGuire (MG) 
energy given by Eq. (10) and a = 1. As expected the 
variational GMG energy is lower than the MG energy 
giving a better determination (upper bound) of the true 
ground state of the many-body system. 

It is important to observe that the axial McGuire wave 
function is invariant by a global translation of the posi- 
tions of the particles. As shown by Calogero and De- 
gasperis [12], and also by Castin and Herzog [14], the 
wave function of Eq. (8) can be normalized by imposing 
that its center of mass Zcm is fixed, namely 



/ 



Siz-Zc^)\f\^ dzi...dzN = N 



The density of particles with respect to Zc 
by 



(14) 

is then given 



S{Zcm)S{zi - dzi ... dZN ■ 



(15) 



Following Calogero and Degasperis [12] we immediately 
find 



7 {-iy+HiN\fc^j>{--/lN\z\/a^) 
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For N ^ 1 this formula can be approximated by 
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(17) 



that is the first term in a 1/A^ expansion of the previ- 
ous expression. Obviously, only for a ~ 1 the Eq. (16) 
gives the ID solitonic profile predicted by Calogero and 
Degasperis [12]. In our variational GMG scheme a is not 
constant and must be determined using the Eq. (11). 



IV. COMPARISON WITH THE 3D GPE 

Here we compare the GMG wave function with the 
fully 3D Hartree wave function, from which one obtains 
the 3D Gross-Pitaevskii equation (GPE). The 3D Hartree 
approximation is obtained by setting 



N 



*(ri,...,rjv)=iVi/2[]V(rO- 



(18) 



Inserting this many-body wave function into Eq. (2) with 
Eq. (4), after integration one finds the Gross-Pitaevskii 
energy functional 



E[^{r)] = N j i)* 



(19) 



By minimizing this energy functional with respect to V'(r) 
with the constraint of the normalization J \'tp\'^d^r = 1 
one finds the familiar mean-field 3D GPE given by 



> V + -{x^ + y^) - 2nj{N - 1)|V'|' 



^ = HiP, (20) 



where /i is the chemical potential fixed by the normal- 
ization. We solve the stationary GPE by using a finite- 
difference Crank-Nicolson predictor-corrector algorithm 
with imaginary time [15]. 
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FIG. 2. Density profile p{z) of the bright soli- 
ton solution obtained with three different theories: 
Gauss-McGuire (GMG), McGuire (MG) and the Gross- 
Pitaevskii equation (GPE). Four values of the interaction 
strength: (a) -/{N - 1) = 0.7; (b) jiN - 1) = 0.9; (c) 
ry{N - 1) = 1; (d) 7(iV - 1) = 1.2. N is the number of 
particles and 7 = 0.001. 

In Fig. 2 we plot the solitonic density profile of the 

GPE and compare it with the profile of GMG theory, i.e. 
Eq. (16) with (7 given by Eq. (11). We insert also the soli- 
tonic profile of the MG theory, i.e. Eq. (16) with a = 1. 
Up to the collapse our GMG solitionic profile is close to 
the GPE one; there are instead relevant deviations with 
respect to the MG soliton (which does not collapse) . Our 
numerical integration of the GPE gives the collapse for 
^{N — 1) ~ 1.35 in agreement with previous computa- 
tions [16] and not too far from the analytical prediction 
of Eq. (13). The GPE critical strength is very close 
to the analytical prediction ^{N — 1) = 4/3 ~ 1.33 of 
the non-polynomial Schrodingcr equation (NPSE) [4,5]. 
However the NPSE is not an exact variational equation: 
it was obtained from the GPE by neglecting some space 
derivatives (for details see [4] and [17]). 

The results of Fig. 2 are obtained by setting 7 = 0.001. 
In Fig. 3 we use the same value of 7, that is relevant for 
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the experiments performed with Bosc-Einstcin conden- 
sates of ^Li atoms [18,19]. Fig. 3 shows the energy per 
particle E/N, the chemical potential fx, the density per 
particle p{0)/N, and also the transverse width a of the 
soliton as a function of the number N of bosons. The 
GMG energy (solid line) is always close to the GPE en- 
ergy (dotted line) and near the collapse the energy differ- 
ence slightly increases. Chemical potentials and densities 
have a similar behavior. Fig. 3 shows that for a large 
number N of bosons the MG theory (dashed lines) dif- 
fers with respect to the other two theories: for N > 750 
deviations are clearly visible. 
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FIG. 3. Bright soliton properties. Comparison among 
different theories: Gauss-McGuire (solid line), McGuire 

(dashed line) and the Gross-Pitaevskii equation (dotted 
line), (a) Energy per particle E/N. (b) Chemical po- 
tential jj,. (c) Axial density per particle p{0)/N at the 
origin, (d) Transverse width a at ^ = 0. A'' is the num- 
ber of particles and 7 = 0.001. 

The panel (d) of Fig. 3 shows that the transverse width 
a decreases by increasing N for both the GPE and the 
GMG theory; however a stays always of the order of one 
up to the collapse. Note that the transverse width a ob- 
tained by using the GPE depends on the axial coordinate 
z. In the plot we choose a at z = 0, that is the lower 
value of the transverse width. 

The variational principle says that the lowest energy 
indicates the most accurate solution. The variational 
principle applies when different variational solutions for 
the same Hamiltonian are compared. In our case, the 
Hamiltonian is given by Eq. (1) with Eq. (4) and the 
three variational solutions are: the GMG one given by 
Eq. (3) and Eq. (8), the MG one given by Eq. (3) with 



(7 = 1 and Eq. (8), and the GPE one given by Eq. (19). 
For N = 1 the three theories give the same value of en- 
ergy: E/N = 1. The top panel of Fig. 3 shows that for 
a large N the GPE gives the lowest energy. We expect 
that for a small N the GMG theory and the MG theory 
are more reliable than the GPE. To verify this predic- 
tion we choose a larger value of 7, namely 7 = 0.01, and 
calculate numerically the GPE energy of Eq. (19) for all 
integer values of N up the the collapse. The results are 
shown in Fig. 4 where we plot the energies for increasing 
values of N. 
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FIG. 4. Energy per particle E/N as a function of the 
number N of bosons. Comparison among different theo- 
ries: Gauss-McGuire (solid line), McGuire (dashed line) 
and the Gross-Pitaevskii equation (filled circles). Here 
7 = 0.01. 

For 1 < A'' < 40 the GMG energy (solid line) and 
the MG energy (dashed line) are indistinguishable but 
both lower than the GPE energy (filled circles). Around 
A^ = 60 the GPE energy becomes lower then the MG en- 
ergy but still higher than the GMG energy. Only around 
A' = 100 the GPE energy becomes smaller than the 
GMG energy. At A = 136 there is the collapse of the 
GPE solitonic solution while the GMG soliton collapses 
at A^ = 153. Just before the collapse of the GPE solu- 
tion, i.e. for A = 135, the relative difference An between 
GPE energy and GMG energy is 0.79%, while Ar be- 
tween GPE energy and MG energy is 2.68% (see also the 
upper panel of Fig. 5). 
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FIG. 5. Energetic relative difference (percentual) 
as a function of the number N of bosons. Circles: Ajj be- 
tween McGuire and Gauss-McGuiro. Filled squares: A/f 
between Gross-Pitaevskii equation and Gauss-McGuire. 
Triangles: A^ between Gross-Pitaevskii equation and 
McGuire. 



V. CONCLUSIONS 

We have introduced a beyond mean-field many-body 
wave function which describes dilute attractive bosons 
under transverse harmonic confinement. This variational 
approach, that we have called Gauss-McGuire theory, 
gives simple analytical formulas for the ground-state en- 
ergy and the solitonic density profile. These formulas 
are in good agreement with the numerical results of the 
3D mean-field theory. By comparing the ground-state 
energies, we have verified that for small numbers of par- 
ticles the Gauss-McGuirc theory is better than the 3D 
mean-field theory while the 3D mean-field theory is more 
reliable for a large number of particles and a small inter- 
atomic strength. 
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Similar results are obtained by setting 7 = 0.001. Here 

the energetic relative difference A^f between GPE and 
MG theory changes sign around N = 270. Instead Aji 
between GPE and GMG theory changes sign around N = 
480. Just before the collapse of the GPE solution, i.e for 
N = 1334, Ar is 2.86% and 1.00% respectively. 

To better study beyond mean-field effects we perform 
numerical calculations for various values of 7 in the inter- 
val [0.01,0.1]. The results are reported in Fig. 5, where 
we plot the energetic relative difference A r for three val- 
ues of the inter-atomic strength 7. The figure shows that 
the energy of the GMG theory is always smaller than 
the energy of the MG theory; in fact, their relative dif- 
ference Afi (curves with circles) is always positive. The 
behavior of the curves in the three panels reveals that by 
growing 7 the range of N reduces and the range of Ar 
increases. The curves with triangles (GPE energy mi- 
nus GMG energy) and also the curves with filled squares 
(GPE energy minus MG energy) show that by fixing the 
number A'^ of particles and increasing the strength 7, the 
beyond mean-field effects become more important. Fig. 
5 shows that for 7 = 0.05 and 7 = 0.1 the GMG theory 
gives a lower energy than the GPE up to the collapse. We 
find that the critical strength 7s beyond which the GMG 
theory is better than the GPE for all iV is 7^ = 0.044; 
in this regime the critical number Nc of particles for the 
collapse predicted by the GMG theory is more reliable 
than the GPE one. 
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